######################################################### #
###Replications: Climate Framing - Lassoplus#######
######################################################### #
#Paper Title: Systematic mapping of climate and environmental framing experiments and re-analysis with computational methods points to omitted interaction bias 
#Authors: Lukas Fesenfeld*,**,�, Liam Beiser-McGrath***, Yixian Sun****, Michael Wicki*, Thomas Bernauer*
#Current affiliations:*University of Bern, Bern, Switzerland; **ETH Z�rich, Z�rich, Switzerland;***Royal Holloway, University of London, London, UK; ****University of Bath, Bath, UK
#�Corresponding author: Lukas Fesenfeld, lukas.fesenfeld@unibe.ch, University of Bern, Fabrikstrasse 8, 3012 Bern, Switzerland


# rm(list = ls())

#install.packages("sparsereg", dependencies = TRUE)
library(tidyverse)
library(foreign)
library(sparsereg)
library(readstata13)
library(fastDummies)
library(dplyr)
library(dotwhisker)
library(imputeTS)
library(broom)
library(estimatr)
library(kableExtra)

###### #
#User written functions#######
##### #

# convert sparsereg output into a nice clean data frame
clean_sparseout <- function(x) {
  
  out_df <- as.data.frame(summary(x))
  out_df$var <- rownames(out_df)
  # keep only if posterior median != 0
  out_df <- out_df %>%
    filter(table.Posterior.Median != 0)
  # classify whether an estimate is a "main" or "interaction" effect
  out_df <- out_df %>%
    mutate(type = ifelse(str_sub(var, 1, 4) == "trea", "Main", "Het"))
  out_df
  
}

# new function for cleaning results

clean_sparsereg2 <- function(out, study = "NA") {
  main_out <- print(out)
  
  out_stats <- as.data.frame(main_out$statistics)
  out_quants <- as.data.frame(main_out$quantiles)
  
  final_out <- cbind(out_stats, out_quants)
  final_out$var <- rownames(final_out)
  
  final_out <- final_out[c("var", "Mean","50%","2.5%","97.5%")]
  final_out$study <- study
  
  final_out
}

# function for creating clean table from sparsereg output

clean_sparsereg_table <- function(x) {
  x %>%
    filter(`50%` != 0) %>% 
    mutate_if(is.numeric, round, 3) %>%
    htmlTable::htmlTable()
}


################################################################################################ #
#ID17: Renewable energy policy design and framing influence public support in the United States#---------
################################################################################################ #
# In Excel file p.2 as ID 13

X_id17 <- as.matrix(data.frame(Republican = id17_lasso$rep, 
                               Democrat = id17_lasso$demo,
                               Independent = id17_lasso$indep, 
                               age = id17_lasso$age, edu = id17_lasso$edu,
                               race = id17_lasso$race, female = id17_lasso$female))

set.seed(12345)
out_id17_lasso <- sparsereg::sparsereg(y = id17_lasso$billsup, 
                                       X = X_id17,
                                       treat = id17_lasso[,paste0("Randomization",1:5,"b")],
                                       scale.type = "TX",
                                       baseline.vec = "Control",
                                       gibbs=1000, burnin=1000)

out_id17_lasso_clean <- clean_sparsereg2(out_id17_lasso, study = "id17")
out_id17_lasso_clean$outcome <- "billsup"

clean_sparsereg_table(out_id17_lasso_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id17_lasso_clean.html")

################################################################################### #
# ID 35: Costs, benefits, and the malleability of public support for “Fracking” ######
################################################################################### #


# recode pid7 (party), in codefile: pid7[1,2,3]= dem5, pid7[5,6,7] = gop5, pid7[4]=ind5
# id35_lasso$republican <- id35_lasso$pid7 ##or what means gop?

id35_lasso <- id35_lasso %>%
  mutate(Republican = ifelse(republican %in% 1:3, 1, 0),
         Democrat = ifelse(republican %in% 5:7, 1, 0),
         Independent = ifelse(republican == 4, 1, 0))

X_id35 <- as.matrix(data.frame(Republican = id35_lasso$Republican, 
                               Democrat = id35_lasso$Democrat,
                               Independent = id35_lasso$Independent, 
                               age = id35_lasso$age, educ = id35_lasso$educ,
                               employ = id35_lasso$employ,
                               white = id35_lasso$white, male = id35_lasso$male))

set.seed(12345)
out_id35_lasso <- sparsereg::sparsereg(y = id35_lasso$DVfrackingsupport, 
                                       X = X_id35,
                                       treat = id35_lasso$treat,
                                       scale.type = "TX",
                                       baseline.vec = "control",
                                       gibbs=1000, burnin=1000)

out_id35_lasso_clean <- clean_sparsereg2(out_id35_lasso, study = "id35")
out_id35_lasso_clean$outcome <- "DVfrackingsupport"

clean_sparsereg_table(out_id35_lasso_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id35_lasso_clean.html")

set.seed(12345)
out_id35_lasso2 <- sparsereg::sparsereg(y = id35_lasso$DVfrackingsupport2, 
                                       X = X_id35,
                                       treat = id35_lasso$treat,
                                       scale.type = "TX",
                                       baseline.vec = "control",
                                       gibbs=1000, burnin=1000)

out_id35_lasso2_clean <- clean_sparsereg2(out_id35_lasso2, study = "id35")
out_id35_lasso2_clean$outcome <- "DVfrackingsupport2"

clean_sparsereg_table(out_id35_lasso2_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id35_lasso2_clean.html")

set.seed(12345)
out_id35_lasso3 <- sparsereg::sparsereg(y = id35_lasso$DVfrackingsupport3, 
                                        X = X_id35,
                                        treat = id35_lasso$treat,
                                        scale.type = "TX",
                                        baseline.vec = "control",
                                        gibbs=1000, burnin=1000)

out_id35_lasso3_clean <- clean_sparsereg2(out_id35_lasso3, study = "id35")
out_id35_lasso3_clean$outcome <- "DVfrackingsupport3"

clean_sparsereg_table(out_id35_lasso3_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id35_lasso3_clean.html")

########################################################################################################################### #
#ID 41: How issue frames shape beliefs about the importance of climate change policy across ideological and partisan groups############
########################################################################################################################### #
# In Excel file p.2 as ID 37

X_id41 <- as.matrix(data.frame(Republican = id41_lasso$party_id_f_Republican,
                               Democrat = id41_lasso$party_id_f_Democrat, 
                               Independent = id41_lasso$party_id_f_Independent,
                               age = id41_lasso$age, female = id41_lasso$female,
                               income = id41_lasso$income, ideology = id41_lasso$ideology))

#Setting seed for replication
set.seed(12345)

#Run LASSOplus for mon-ranked climate importance 

import_clim_out_id41_lasso <- sparsereg::sparsereg(y = id41_lasso$importance_climate, 
                                                   X = X_id41,
                                                   treat = id41_lasso$treat,
                                                   scale.type = "TX",
                                                   baseline.vec = "control",
                                                   gibbs=1000, burnin=1000)

out_id41_lasso_clean_imp <- clean_sparsereg2(import_clim_out_id41_lasso, study = "id41")
out_id41_lasso_clean_imp$outcome <- "importance_climate"

clean_sparsereg_table(out_id41_lasso_clean_imp %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id41_lasso_clean_imp.html")

#Run LASSOplus for ranked climate importance 


ranked_import_clim_out_id41_lasso <- sparsereg::sparsereg(y = id41_lasso$rank_climate_rev, 
                                                          X = X_id41,
                                                          treat = id41_lasso$treat,
                                                          scale.type = "TX",
                                                          baseline.vec = "control",
                                                          gibbs=1000, burnin=1000)

out_id41_lasso_clean_imprank <- clean_sparsereg2(ranked_import_clim_out_id41_lasso, study = "id41")
out_id41_lasso_clean_imprank$outcome <- "rank_climate_rev"

clean_sparsereg_table(out_id41_lasso_clean_imprank %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id41_lasso_clean_imprank.html")



####################################################################################################################################### #
# ID 57: Does the label really matter? Evidence that the US public continues to doubt "global warming" more than "climate change" ######
####################################################################################################################################### #

# comparison of effect of cc vs. gw framing by PID (binary outcome: ccgwhappeningbin)

X_id57 <- as.matrix(data.frame(Republican = id57_lasso$republicanbin,
                               Democrat = id57_lasso$democratnbin, 
                               Independent = id57_lasso$politmiddlebin,
                               age = id57_lasso$age, female = id57_lasso$female,
                               educ = id57_lasso$educ))

out_id57_lasso <- sparsereg::sparsereg(y = id57_lasso$ccgwhappeningbin, 
                                       X = X_id57,
                                       treat = id57_lasso$treat,
                                       scale.type = "TX",
                                       baseline.vec = "globalwarming",
                                       gibbs=1000, burnin=1000)

out_id57_lasso_clean <- clean_sparsereg2(out_id57_lasso, study = "id57")
out_id57_lasso_clean$outcome <- "ccgwhappeningbin"

clean_sparsereg_table(out_id57_lasso_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id57_lasso_clean.html")

########################################################################################## #
# ID 71: Doing What Others Do: Norms, Science, and Collective Action on Global Warming ######
########################################################################################## #

# experiment 1
names(id71_10_lasso)

# generate party dummies
# negative is republican according to codebook
# no other control vars

id71_10_lasso <- id71_10_lasso %>%
  mutate(pidlab = case_when(partyid < -0.4 ~ "Republican",
                            partyid >= -0.4 & partyid < 0.4 ~ "Independent",
                            partyid >= 0.4 ~ "Democrat")) %>%
  fastDummies::dummy_cols(select_columns = "pidlab")


X_id71_10 <- as.matrix(data.frame(Republican = id71_10_lasso$pidlab_Republican,
                                  Democrat = id71_10_lasso$pidlab_Democrat, 
                                  Independent = id71_10_lasso$pidlab_Independent))

# outcomes: ce_ppi, ce_er, global_warming, human_induced, support_cap, smallercar

set.seed(12345)
out_id71_10_lasso_ce_ppi <- sparsereg::sparsereg(y = id71_10_lasso$ce_ppi, 
                                       X = X_id71_10,
                                       treat = id71_10_lasso$treat,
                                       scale.type = "TX",
                                       baseline.vec = "control",
                                       gibbs=1000, burnin=1000)

out_id71_10_lasso_ce_ppi_clean <- clean_sparsereg2(out_id71_10_lasso_ce_ppi, study = "id71")
out_id71_10_lasso_ce_ppi_clean$outcome <- "ce_ppi"



clean_sparsereg_table(out_id71_10_lasso_ce_ppi_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_10_lasso_ce_ppi_clean.html")

set.seed(12345)
out_id71_10_lasso_ce_er <- sparsereg::sparsereg(y = id71_10_lasso$ce_er, 
                                                 X = X_id71_10,
                                                 treat = id71_10_lasso$treat,
                                                 scale.type = "TX",
                                                 baseline.vec = "control",
                                                 gibbs=1000, burnin=1000)
out_id71_10_lasso_ce_er_clean <- clean_sparsereg2(out_id71_10_lasso_ce_er, study = "id71")
out_id71_10_lasso_ce_er_clean$outcome <- "ce_er"

clean_sparsereg_table(out_id71_10_lasso_ce_er_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_10_lasso_ce_er_clean.html")

set.seed(12345)
out_id71_10_lasso_global_warming <- sparsereg::sparsereg(y = id71_10_lasso$global_warming, 
                                                 X = X_id71_10,
                                                 treat = id71_10_lasso$treat,
                                                 scale.type = "TX",
                                                 baseline.vec = "control",
                                                 gibbs=1000, burnin=1000)
out_id71_10_lasso_global_warming_clean <- clean_sparsereg2(out_id71_10_lasso_global_warming, study = "id71")
out_id71_10_lasso_global_warming_clean$outcome <- "global_warming"

clean_sparsereg_table(out_id71_10_lasso_global_warming_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_10_lasso_global_warming_clean.html")


set.seed(12345)
out_id71_10_lasso_human_induced <- sparsereg::sparsereg(y = id71_10_lasso$human_induced, 
                                                 X = X_id71_10,
                                                 treat = id71_10_lasso$treat,
                                                 scale.type = "TX",
                                                 baseline.vec = "control",
                                                 gibbs=1000, burnin=1000)
out_id71_10_lasso_human_induced_clean <- clean_sparsereg2(out_id71_10_lasso_human_induced, study = "id71")
out_id71_10_lasso_human_induced_clean$outcome <- "human_induced"

clean_sparsereg_table(out_id71_10_lasso_human_induced_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_10_lasso_human_induced_clean.html")

set.seed(12345)
out_id71_10_lasso_support_cap <- sparsereg::sparsereg(y = id71_10_lasso$support_cap, 
                                                 X = X_id71_10,
                                                 treat = id71_10_lasso$treat,
                                                 scale.type = "TX",
                                                 baseline.vec = "control",
                                                 gibbs=1000, burnin=1000)
out_id71_10_lasso_support_cap_clean <- clean_sparsereg2(out_id71_10_lasso_support_cap, study = "id71")
out_id71_10_lasso_support_cap_clean$outcome <- "support_cap"

clean_sparsereg_table(out_id71_10_lasso_support_cap_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_10_lasso_support_cap_clean.html")


set.seed(12345)
out_id71_10_lasso_smallercar <- sparsereg::sparsereg(y = id71_10_lasso$smallercar, 
                                                 X = X_id71_10,
                                                 treat = id71_10_lasso$treat,
                                                 scale.type = "TX",
                                                 baseline.vec = "control",
                                                 gibbs=1000, burnin=1000)
out_id71_10_lasso_smallercar_clean <- clean_sparsereg2(out_id71_10_lasso_smallercar, study = "id71")
out_id71_10_lasso_smallercar_clean$outcome <- "smallercar"

clean_sparsereg_table(out_id71_10_lasso_smallercar_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_10_lasso_smallercar_clean.html")


# experiment 2
names(id71_11_lasso)



# generate party dummies
# negative is republican according to codebook

id71_11_lasso <- id71_11_lasso %>%
  mutate(pidlab = case_when(democrat < -0.4 ~ "Republican",
                            democrat >= -0.4 & democrat < 0.4 ~ "Independent",
                            democrat >= 0.4 ~ "Democrat")) %>%
  fastDummies::dummy_cols(select_columns = "pidlab")


X_id71_11 <- as.matrix(data.frame(Republican = id71_11_lasso$pidlab_Republican,
                               Democrat = id71_11_lasso$pidlab_Democrat, 
                               Independent = id71_11_lasso$pidlab_Independent,
                               age = id71_11_lasso$age, female = id71_11_lasso$female,
                               liberal = id71_11_lasso$liberal))

# outcomes: GWBelief, GWHuman, EmissionCap, PersonalAction

set.seed(12345)
out_id71_11_lasso_GWBelief <- sparsereg::sparsereg(y = id71_11_lasso$GWBelief, 
                                                   X = X_id71_11,
                                                   treat = id71_11_lasso$treat,
                                                   scale.type = "TX",
                                                   baseline.vec = "control",
                                                   gibbs=1000, burnin=1000)
out_id71_11_lasso_GWBelief_clean <- clean_sparsereg2(out_id71_11_lasso_GWBelief, study = "id71")
out_id71_11_lasso_GWBelief_clean$outcome <- "GWBelief"

clean_sparsereg_table(out_id71_11_lasso_GWBelief_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_11_lasso_GWBelief_clean.html")

set.seed(12345)
out_id71_11_lasso_GWHuman <- sparsereg::sparsereg(y = id71_11_lasso$GWHuman, 
                                                  X = X_id71_11,
                                                  treat = id71_11_lasso$treat,
                                                  scale.type = "TX",
                                                  baseline.vec = "control",
                                                  gibbs=1000, burnin=1000)
out_id71_11_lasso_GWHuman_clean <- clean_sparsereg2(out_id71_11_lasso_GWHuman, study = "id71")
out_id71_11_lasso_GWHuman_clean$outcome <- "GWHuman"

clean_sparsereg_table(out_id71_11_lasso_GWHuman_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_11_lasso_GWHuman_clean.html")


set.seed(12345)
out_id71_11_lasso_EmissionCap <- sparsereg::sparsereg(y = id71_11_lasso$EmissionCap, 
                                                      X = X_id71_11,
                                                      treat = id71_11_lasso$treat,
                                                      scale.type = "TX",
                                                      baseline.vec = "control",
                                                      gibbs=1000, burnin=1000)
out_id71_11_lasso_EmissionCap_clean <- clean_sparsereg2(out_id71_11_lasso_EmissionCap, study = "id71")
out_id71_11_lasso_EmissionCap_clean$outcome <- "EmissionCap"

clean_sparsereg_table(out_id71_11_lasso_EmissionCap_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_11_lasso_EmissionCap_clean.html")


set.seed(12345)
out_id71_11_lasso_PersonalAction <- sparsereg::sparsereg(y = id71_11_lasso$PersonalAction, 
                                                         X = X_id71_11,
                                                         treat = id71_11_lasso$treat,
                                                         scale.type = "TX",
                                                         baseline.vec = "control",
                                                         gibbs=1000, burnin=1000)
out_id71_11_lasso_PersonalAction_clean <- clean_sparsereg2(out_id71_11_lasso_PersonalAction, study = "id71")
out_id71_11_lasso_PersonalAction_clean$outcome <- "PersonalAction"

clean_sparsereg_table(out_id71_11_lasso_PersonalAction_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id71_11_lasso_PersonalAction_clean.html")

########################################################################################################################################### #
# ID 73: The impact of elite frames and motivated reasoning on beliefs in a global warming conspiracy: The promise and limits of trust ######
########################################################################################################################################### #

# note id73_12_lasso and id73_13_lasso should be the same (difference in Stata format)
# so will just use id73_13_lasso

id73_13_lasso <- id73_13_lasso %>%
  mutate(Democrat = ifelse(republican %in% 1:3, 1, 0),
         Republican = ifelse(republican %in% 5:7, 1, 0),
         Independent = ifelse(republican == 4, 1, 0))

X_id73_13 <- as.matrix(data.frame(Republican = id73_13_lasso$Republican,
                                  Democrat = id73_13_lasso$Democrat, 
                                  Independent = id73_13_lasso$Independent,
                                  age = id73_13_lasso$age, female = id73_13_lasso$female,
                                  income = id73_13_lasso$income,
                                  white = id73_13_lasso$white,
                                  educ_alt = id73_13_lasso$educ_alt))

set.seed(12345)
out_id73_13_lasso <- sparsereg::sparsereg(y = id73_13_lasso$ccgwishoax,
                                          X = X_id73_13,
                                          treat = id73_13_lasso$treat,
                                          scale.type = "TX",
                                          baseline.vec = "GlobalWarming",
                                          gibbs=1000, burnin=1000)
out_id73_13_lasso_clean <- clean_sparsereg2(out_id73_13_lasso, study = "id73")
out_id73_13_lasso_clean$outcome <- "ccgwishoax"

clean_sparsereg_table(out_id73_13_lasso_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id73_13_lasso_clean.html")


####################################################################################################### #
# ID 74: A Dirty Word or a Dirty World? Attribute Framing, Political Affiliation, and Query Theory ######
####################################################################################################### #

id74_1_lasso <- id74_1_lasso %>%
  mutate(Democrat = ifelse(republican1 == -1, 1, 0),
         Independent = ifelse(republican1 == 0, 1, 0),
         Republican = ifelse(republican1 == 1, 1, 0))

X_id74_1 <- as.matrix(data.frame(Democrat = id74_1_lasso$Democrat,
                                       Independent = id74_1_lasso$Independent,
                                       Republican = id74_1_lasso$Republican,
                                       AGE = id74_1_lasso$AGE,
                                       education = id74_1_lasso$education,
                                       income = id74_1_lasso$income,
                                       female = id74_1_lasso$M0F1))

# outcome: co_avg, ps_avg, mm_avg, cr_avg

set.seed(12345)
out_id74_1_lasso_co_avg <- sparsereg::sparsereg(y = id74_1_lasso$co_avg,
                                          X = X_id74_1,
                                          treat = id74_1_lasso$treat,
                                          scale.type = "TX",
                                          baseline.vec = "offset",
                                          gibbs=1000, burnin=1000)
out_id74_1_lasso_co_avg_clean <- clean_sparsereg2(out_id74_1_lasso_co_avg, study = "id74")
out_id74_1_lasso_co_avg_clean$outcome <- "co_avg"

clean_sparsereg_table(out_id74_1_lasso_co_avg_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id74_1_lasso_co_avg_clean.html")


set.seed(12345)
out_id74_1_lasso_ps_avg <- sparsereg::sparsereg(y = id74_1_lasso$ps_avg,
                                                X = X_id74_1,
                                                treat = id74_1_lasso$treat,
                                                scale.type = "TX",
                                                baseline.vec = "offset",
                                                gibbs=1000, burnin=1000)
out_id74_1_lasso_ps_avg_clean <- clean_sparsereg2(out_id74_1_lasso_ps_avg, study = "id74")
out_id74_1_lasso_ps_avg_clean$outcome <- "ps_avg"

clean_sparsereg_table(out_id74_1_lasso_ps_avg_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id74_1_lasso_ps_avg_clean.html")


set.seed(12345)
out_id74_1_lasso_mm_avg <- sparsereg::sparsereg(y = id74_1_lasso$mm_avg,
                                                X = X_id74_1,
                                                treat = id74_1_lasso$treat,
                                                scale.type = "TX",
                                                baseline.vec = "offset",
                                                gibbs=1000, burnin=1000)
out_id74_1_lasso_mm_avg_clean <- clean_sparsereg2(out_id74_1_lasso_mm_avg, study = "id74")
out_id74_1_lasso_mm_avg_clean$outcome <- "mm_avg"

clean_sparsereg_table(out_id74_1_lasso_mm_avg_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id74_1_lasso_mm_avg_clean.html")


set.seed(12345)
out_id74_1_lasso_cr_avg <- sparsereg::sparsereg(y = id74_1_lasso$cr_avg,
                                                X = X_id74_1,
                                                treat = id74_1_lasso$treat,
                                                scale.type = "TX",
                                                baseline.vec = "offset",
                                                gibbs=1000, burnin=1000)
out_id74_1_lasso_cr_avg_clean <- clean_sparsereg2(out_id74_1_lasso_cr_avg, study = "id74")
out_id74_1_lasso_cr_avg_clean$outcome <- "cr_avg"

clean_sparsereg_table(out_id74_1_lasso_cr_avg_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id74_1_lasso_cr_avg_clean.html")


# experiment 2
# outcome: choice, pref, mand

X_id74_2 <- as.matrix(data.frame(Democrat = id74_2_lasso$Politicalaffiliation_Democrat,
                                       Independent = id74_2_lasso$Politicalaffiliation_Independent,
                                       Republican = id74_2_lasso$Politicalaffiliation_Republican,
                                       HighestEducation = id74_2_lasso$HighestEducation,
                                       age = id74_2_lasso$age, white = id74_2_lasso$white,
                                       Conservatism = id74_2_lasso$Conservatism))

set.seed(12345)
out_id74_2_lasso_choice <- sparsereg::sparsereg(y = id74_2_lasso$choice,
                                                X = X_id74_2,
                                                treat = id74_2_lasso$treat,
                                                scale.type = "TX",
                                                baseline.vec = "offset",
                                                gibbs=1000, burnin=1000)
out_id74_2_lasso_choice_clean <- clean_sparsereg2(out_id74_2_lasso_choice, study = "id74")
out_id74_2_lasso_choice_clean$outcome <- "choice"

clean_sparsereg_table(out_id74_2_lasso_choice_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id74_2_lasso_choice_clean.html")


set.seed(12345)
out_id74_2_lasso_pref <- sparsereg::sparsereg(y = id74_2_lasso$pref,
                                                X = X_id74_2,
                                                treat = id74_2_lasso$treat,
                                                scale.type = "TX",
                                                baseline.vec = "offset",
                                                gibbs=1000, burnin=1000)
out_id74_2_lasso_pref_clean <- clean_sparsereg2(out_id74_2_lasso_pref, study = "id74")
out_id74_2_lasso_pref_clean$outcome <- "pref"

clean_sparsereg_table(out_id74_2_lasso_pref_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id74_2_lasso_pref_clean.html")


set.seed(12345)
out_id74_2_lasso_mand <- sparsereg::sparsereg(y = id74_2_lasso$mand,
                                                X = X_id74_2,
                                                treat = id74_2_lasso$treat,
                                                scale.type = "TX",
                                                baseline.vec = "offset",
                                                gibbs=1000, burnin=1000)
out_id74_2_lasso_mand_clean <- clean_sparsereg2(out_id74_2_lasso_mand, study = "id74")
out_id74_2_lasso_mand_clean$outcome <- "mand"

clean_sparsereg_table(out_id74_2_lasso_mand_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id74_2_lasso_mand_clean.html")


####################################################################################################################### #
#ID 83: Do partisanship and politicization undermine the impact of a scientific consensus message about  climate change?#####
####################################################################################################################### ##

X_id83 <- as.matrix(data.frame(Democrat = id83_lasso_part$dem, 
                               Republican = id83_lasso_part$rep, knowdummy = id83_lasso_part$knowdummy,
                               inform = id83_lasso_part$inform, politi = id83_lasso_part$politi,
                               warning = id83_lasso_part$warning, correction = id83_lasso_part$correction,
                               female = id83_lasso_part$female, agerange = id83_lasso_part$agerange,
                               education = id83_lasso_part$education, income = id83_lasso_part$income,
                               minority = id83_lasso_part$minority))

set.seed(12345)

out_id83_lasso_part_scienceagree_main <- sparsereg::sparsereg(y = id83_lasso_part$scienceagree, 
                                                              X = X_id83,
                                                              treat = id83_lasso_part$treatment_f,
                                                              scale.type = "TX",
                                                              baseline.vec = "Control",
                                                              gibbs=1000, burnin=1000)

out_id83_lasso_part_scienceagree_main_clean <- clean_sparsereg2(out_id83_lasso_part_scienceagree_main, study = "id83")
out_id83_lasso_part_scienceagree_main_clean$outcome <- "scienceagree"

clean_sparsereg_table(out_id83_lasso_part_scienceagree_main_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id83_lasso_part_scienceagree_main_clean.html")

set.seed(12345)

out_id83_lasso_part_policy_main <- sparsereg::sparsereg(y = id83_lasso_part$policy, 
                                                        X = X_id83,
                                                        treat = id83_lasso_part$treatment_f,
                                                        scale.type = "TX",
                                                        baseline.vec = "Control",
                                                        gibbs=1000, burnin=1000)

out_id83_lasso_part_policy_main_clean <- clean_sparsereg2(out_id83_lasso_part_policy_main, study = "id83")
out_id83_lasso_part_policy_main_clean$outcome <- "policy"

clean_sparsereg_table(out_id83_lasso_part_policy_main_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id83_lasso_part_policy_main_clean.html")


set.seed(12345)

out_id83_lasso_part_ccexphuman_main <- sparsereg::sparsereg(y = id83_lasso_part$ccexphuman, 
                                                            X = X_id83,
                                                            treat = id83_lasso_part$treatment_f,
                                                            scale.type = "TX",
                                                            baseline.vec = "Control",
                                                            gibbs=1000, burnin=1000)

out_id83_lasso_part_ccexphuman_main_clean <- clean_sparsereg2(out_id83_lasso_part_ccexphuman_main, study = "id83")
out_id83_lasso_part_ccexphuman_main_clean$outcome <- "ccexphuman"

clean_sparsereg_table(out_id83_lasso_part_ccexphuman_main_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id83_lasso_part_ccexphuman_main_clean.html")


########################################################################################################### #
# ID 86: "Global warming" or "climate change"?: Whether the planet is warming depends on question wording ######
########################################################################################################### #


X_id86 <- as.matrix(data.frame(Democrat = id86_lasso$p_Democrat, 
                               Republican = id86_lasso$p_Republican,
                               Independent = id86_lasso$p_Independent,
                               female = id86_lasso$female, age = id86_lasso$age,
                               educ = id86_lasso$educ, white = id86_lasso$white,
                               income = id86_lasso$familyincome))

set.seed(12345)
out_id86_lasso <- sparsereg::sparsereg(y = id86_lasso$gwcchappen, 
                                       X = X_id86,
                                       treat = id86_lasso$treat,
                                       scale.type = "TX",
                                       baseline.vec = "globalwarming",
                                       gibbs=1000, burnin=1000)
out_id86_lasso_clean <- clean_sparsereg2(out_id86_lasso, study = "id86")
out_id86_lasso_clean$outcome <- "gwcchappen"

clean_sparsereg_table(out_id86_lasso_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id86_lasso_clean.html")

 

############################################################################################################################################# #
#ID 109: Credibility, communication, and climate change: How lifestyle inconsistency and do-gooder derogation impact decarbonization advocacy#
############################################################################################################################################# #

# not clear if really subgroup effects here
# X_id109 <- as.matrix(data.frame(Pol = id109_lasso$Pol, 
#                                 Income = id109_lasso$Income,
#                                 Edu = id109_lasso$Edu,
#                                 female = id109_lasso$female, Age = id109_lasso$Age))
# 



####################################################################################################################### #
# ID 115: Economic losses or environmental gains? Framing effects on public support for environmental management ######
####################################################################################################################### #

id115_lasso <- id115_lasso %>%
  mutate(pid_label = case_when(partyid == 1 ~ "Republican",
                               partyid == 2 ~ "Democrat",
                               partyid == 3 ~ "Independent",
                               partyid == 4 ~ "Other")) %>%
  fastDummies::dummy_cols(select_columns = "pid_label")

X_id115 <- as.matrix(data.frame(Republican = id115_lasso$pid_label_Republican,
                                Democrat = id115_lasso$pid_label_Democrat,
                                Independent = id115_lasso$pid_label_Independent,
                                female = id115_lasso$female, age = id115_lasso$age,
                                white = id115_lasso$white, black = id115_lasso$black,
                                amindian = id115_lasso$am_indian, asian = id115_lasso$asian,
                                pacislander = id115_lasso$pac_islander, latino = id115_lasso$latino,
                                raceother = id115_lasso$race_other, hhincome = id115_lasso$hhincome))

set.seed(12345)
out_id115_lasso <- sparsereg::sparsereg(y = id115_lasso$support, 
                                       X = X_id115,
                                       treat = id115_lasso$treat,
                                       scale.type = "TX",
                                       baseline.vec = "control",
                                       gibbs=1000, burnin=1000)

out_id115_lasso_clean <- clean_sparsereg2(out_id115_lasso, study = "id115")
out_id115_lasso_clean$outcome <- "support"

clean_sparsereg_table(out_id115_lasso_clean %>%
                        dplyr::select(-var, -study)) %>%
  save_kable("./01_Data Analysis/output/out_id115_lasso_clean.html")



######################## 3
#### OLS Analysis
########################## #



# ID17 - OLS Sub-group
out_id17_ols <- id17_lasso %>%
  mutate(pid = case_when(rep == 1 ~ "Republican",
                         demo == 1 ~ "Democrat",
                         indep == 1 ~ "Independent")) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("billsup ~ ",
                                   paste0("Randomization",1:5,"b", collapse = " + "))),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# ID35 - OLS Sub-group
out_id35_ols <- id35_lasso %>%
  mutate(pid = case_when(republican %in% 1:3 ~ "Republican",
                         republican == 4 ~ "Independent",
                         republican %in% 5:7 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("DVfrackingsupport ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id35_ols2 <- id35_lasso %>%
  mutate(pid = case_when(republican %in% 1:3 ~ "Republican",
                         republican == 4 ~ "Independent",
                         republican %in% 5:7 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("DVfrackingsupport2 ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id35_ols3 <- id35_lasso %>%
  mutate(pid = case_when(republican %in% 1:3 ~ "Republican",
                         republican == 4 ~ "Independent",
                         republican %in% 5:7 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("DVfrackingsupport3 ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# ID41 - OLS Sub-group
out_id41_ols_imp <- id41_lasso %>%
  mutate(pid = case_when(party_id_f_Republican == 1 ~ "Republican",
                         party_id_f_Independent == 1 ~ "Independent",
                         party_id_f_Democrat == 1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("importance_climate ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id41_ols_imprank <- id41_lasso %>%
  mutate(pid = case_when(party_id_f_Republican == 1 ~ "Republican",
                         party_id_f_Independent == 1 ~ "Independent",
                         party_id_f_Democrat == 1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("rank_climate_rev ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# ID57 - OLS Sub-group
out_id57_ols <- id57_lasso %>%
  mutate(pid = case_when(republicanbin == 1 ~ "Republican",
                         politmiddlebin == 1 ~ "Independent",
                         democratnbin == 1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "globalwarming", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("ccgwhappeningbin ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# ID71 - OLS Sub-group
out_id71_10_ols_ce_ppi <- id71_10_lasso %>%
  mutate(pid = case_when(partyid < -0.4 ~ "Republican",
                         partyid >= -0.4 & partyid < 0.4 ~ "Independent",
                         partyid >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("ce_ppi ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id71_10_ols_ce_er <- id71_10_lasso %>%
  mutate(pid = case_when(partyid < -0.4 ~ "Republican",
                         partyid >= -0.4 & partyid < 0.4 ~ "Independent",
                         partyid >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("ce_er ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id71_10_ols_global_warming <- id71_10_lasso %>%
  mutate(pid = case_when(partyid < -0.4 ~ "Republican",
                         partyid >= -0.4 & partyid < 0.4 ~ "Independent",
                         partyid >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("global_warming ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id71_10_ols_human_induced <- id71_10_lasso %>%
  mutate(pid = case_when(partyid < -0.4 ~ "Republican",
                         partyid >= -0.4 & partyid < 0.4 ~ "Independent",
                         partyid >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("human_induced ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id71_10_ols_support_cap <- id71_10_lasso %>%
  mutate(pid = case_when(partyid < -0.4 ~ "Republican",
                         partyid >= -0.4 & partyid < 0.4 ~ "Independent",
                         partyid >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("support_cap ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id71_10_ols_smallercar <- id71_10_lasso %>%
  mutate(pid = case_when(partyid < -0.4 ~ "Republican",
                         partyid >= -0.4 & partyid < 0.4 ~ "Independent",
                         partyid >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("smallercar ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# experiment 2

out_id71_11_ols_GWBelief <- id71_11_lasso %>%
  mutate(pid = case_when(democrat < -0.4 ~ "Republican",
                         democrat >= -0.4 & democrat < 0.4 ~ "Independent",
                         democrat >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("GWBelief ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id71_11_ols_GWHuman <- id71_11_lasso %>%
  mutate(pid = case_when(democrat < -0.4 ~ "Republican",
                         democrat >= -0.4 & democrat < 0.4 ~ "Independent",
                         democrat >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("GWHuman ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id71_11_ols_EmissionCap <- id71_11_lasso %>%
  mutate(pid = case_when(democrat < -0.4 ~ "Republican",
                         democrat >= -0.4 & democrat < 0.4 ~ "Independent",
                         democrat >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("EmissionCap ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id71_11_ols_PersonalAction <- id71_11_lasso %>%
  mutate(pid = case_when(democrat < -0.4 ~ "Republican",
                         democrat >= -0.4 & democrat < 0.4 ~ "Independent",
                         democrat >= 0.4 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("PersonalAction ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# ID73 - OLS Sub-group
out_id73_13_ols <- id73_13_lasso %>%
  mutate(pid = case_when(republican %in% 1:3 ~ "Democrat",
                         republican == 4 ~ "Independent",
                         republican %in% 5:7 ~ "Republican")) %>%
  mutate(treat = fct_relevel(treat, "GlobalWarming", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("ccgwishoax ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# ID74 - OLS Sub-group
out_id74_1_ols_co_avg <- id74_1_lasso %>%
  mutate(pid = case_when(republican1 %in% 1 ~ "Republican",
                         republican1 == 0 ~ "Independent",
                         republican1 %in% -1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "offset", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("co_avg ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id74_1_ols_ps_avg <- id74_1_lasso %>%
  mutate(pid = case_when(republican1 %in% 1 ~ "Republican",
                         republican1 == 0 ~ "Independent",
                         republican1 %in% -1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "offset", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("ps_avg ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id74_1_ols_mm_avg <- id74_1_lasso %>%
  mutate(pid = case_when(republican1 %in% 1 ~ "Republican",
                         republican1 == 0 ~ "Independent",
                         republican1 %in% -1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "offset", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("mm_avg ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id74_1_ols_cr_avg <- id74_1_lasso %>%
  mutate(pid = case_when(republican1 %in% 1 ~ "Republican",
                         republican1 == 0 ~ "Independent",
                         republican1 %in% -1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "offset", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("cr_avg ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# experiment 2

out_id74_2_ols_choice <- id74_2_lasso %>%
  mutate(pid = case_when(Politicalaffiliation_Republican == 1 ~ "Republican",
                         Politicalaffiliation_Independent == 1 ~ "Independent",
                         Politicalaffiliation_Democrat == 1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "offset", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("choice ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id74_2_ols_pref <- id74_2_lasso %>%
  mutate(pid = case_when(Politicalaffiliation_Republican == 1 ~ "Republican",
                         Politicalaffiliation_Independent == 1 ~ "Independent",
                         Politicalaffiliation_Democrat == 1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "offset", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("pref ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))


out_id74_2_ols_mand <- id74_2_lasso %>%
  mutate(pid = case_when(Politicalaffiliation_Republican == 1 ~ "Republican",
                         Politicalaffiliation_Independent == 1 ~ "Independent",
                         Politicalaffiliation_Democrat == 1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "offset", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("mand ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# ID83 - OLS Sub-group
out_id83_ols_scienceagree <- id83_lasso %>%
  mutate(pid = case_when(rep == 1 ~ "Republican",
                         dem == 1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treatment_f, "Control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("scienceagree ~ ",
                                   "treatment_f")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id83_ols_policy <- id83_lasso %>%
  mutate(pid = case_when(rep == 1 ~ "Republican",
                         dem == 1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treatment_f, "Control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("policy ~ ",
                                   "treatment_f")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

out_id83_ols_ccexphuman <- id83_lasso %>%
  mutate(pid = case_when(rep == 1 ~ "Republican",
                         dem == 1 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treatment_f, "Control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("ccexphuman ~ ",
                                   "treatment_f")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

# ID86 - OLS Sub-group
out_id86_ols <- id86_lasso %>%
  mutate(pid = case_when(p_Republican == 1  ~ "Republican",
                         p_Independent == 1  ~ "Independent",
                         p_Democrat == 1  ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "globalwarming", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("gwcchappen ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))


# ID115 - OLS Sub-group
out_id115_ols <- id115_lasso %>%
  mutate(pid = case_when(partyid == 1 ~ "Republican",
                         partyid == 3 ~ "Independent",
                         partyid == 2 ~ "Democrat")) %>%
  mutate(treat = fct_relevel(treat, "control", after = 0)) %>%
  group_by(pid) %>%
  do(tidy(lm_robust(formula(paste0("support ~ ",
                                   "treat")),
                    data = .))) %>%
  filter(term != "(Intercept)") %>%
  filter(!is.na(pid))

#### cleaning up results

ols_comb_out <- bind_rows("id115" = out_id115_ols,                       
                          "id17" = out_id17_ols,                        
                          "id35" = out_id35_ols,                        
                          "id35" = out_id35_ols2,                       
                          "id35" = out_id35_ols3,                       
                          "id41" = out_id41_ols_imp,                   
                          "id41" = out_id41_ols_imprank,               
                          "id57" = out_id57_ols,                        
                          "id71" = out_id71_10_ols_ce_er,               
                          "id71" = out_id71_10_ols_ce_ppi,              
                          "id71" = out_id71_10_ols_global_warming,      
                          "id71" = out_id71_10_ols_human_induced,       
                          "id71" = out_id71_10_ols_smallercar,          
                          "id71" = out_id71_10_ols_support_cap,         
                          "id71" = out_id71_11_ols_EmissionCap,         
                          "id71" = out_id71_11_ols_GWBelief,            
                          "id71" = out_id71_11_ols_GWHuman,             
                          "id71" = out_id71_11_ols_PersonalAction,      
                          "id73" = out_id73_13_ols,                     
                          "id74" = out_id74_1_ols_co_avg,               
                          "id74" = out_id74_2_ols_choice,               
                          "id74" = out_id74_2_ols_mand,                 
                          "id74" = out_id74_2_ols_pref,                 
                          "id83" = out_id83_ols_ccexphuman,   
                          "id83" = out_id83_ols_policy,       
                          "id83" = out_id83_ols_scienceagree, 
                          "id86" = out_id86_ols,
                          .id = "study")

ols_comb_out <- ols_comb_out %>%
  mutate(study_lab = case_when(study == "id17" ~ "Stokes and Warshaw (2017)",
                               study == "id35" ~ "Christenson and Goldfarb (2017)",
                               study == "id41" ~ "Singh and Swanson (2017)",
                               study == "id57" ~ "Schuldt and Enns (2017)",
                               study == "id71" ~ "Bolsen and Leeper (2014)",
                               study == "id73" ~ "Saunders (2017)",
                               study == "id74" ~ "Hardisty and Johnson (2010)",
                               study == "id83" ~ "Bolsen and Druckmann (2018)",
                               study == "id86" ~ "Schuldt and Konrath (2011)",
                               study == "id115" ~ "DeGolia and Hiroyasu (2019)"))


ols_comb_out <- ols_comb_out %>%
  mutate(treatment_lab = term) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "treatment_f", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "treatment", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "treat", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization1b", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization2b", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization3b", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization4b", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization5b", ""))

ols_comb_out <- ols_comb_out %>%
  dplyr::rename("Moderator" = "pid")

####### Figures of effect heterogeneity by PID #######

# combine results together
comb_out_clean <- bind_rows(out_id115_lasso_clean,                       
                            out_id17_lasso_clean,                        
                            out_id35_lasso_clean,                        
                            out_id35_lasso2_clean,                       
                            out_id35_lasso3_clean,                       
                            out_id41_lasso_clean_imp,                   
                            out_id41_lasso_clean_imprank,               
                            out_id57_lasso_clean,                        
                            out_id71_10_lasso_ce_er_clean,               
                            out_id71_10_lasso_ce_ppi_clean,              
                            out_id71_10_lasso_global_warming_clean,      
                            out_id71_10_lasso_human_induced_clean,       
                            out_id71_10_lasso_smallercar_clean,          
                            out_id71_10_lasso_support_cap_clean,         
                            out_id71_11_lasso_EmissionCap_clean,         
                            out_id71_11_lasso_GWBelief_clean,            
                            out_id71_11_lasso_GWHuman_clean,             
                            out_id71_11_lasso_PersonalAction_clean,      
                            out_id73_13_lasso_clean,                     
                            out_id74_1_lasso_co_avg_clean,               
                            out_id74_2_lasso_choice_clean,               
                            out_id74_2_lasso_mand_clean,                 
                            out_id74_2_lasso_pref_clean,                 
                            out_id83_lasso_part_ccexphuman_main_clean,   
                            out_id83_lasso_part_policy_main_clean,       
                            out_id83_lasso_part_scienceagree_main_clean, 
                            out_id86_lasso_clean)

# saveRDS(comb_out_clean, file = "./03_Replication_Data_Files/comb_out_clean.RDS")
# saveRDS(comb_out_clean, file = "./03_Replication_Data_Files/comb_out_clean210304.RDS")
# comb_out_clean <- readRDS(file = "./03_Replication_Data_Files/comb_out_clean.RDS")
comb_out_clean <- readRDS(file = "./03_Replication_Data_Files/comb_out_clean210304.RDS")


set.seed(1)

# republican moderator
rephet_results <- comb_out_clean %>%
  filter(grepl("Republican x", var))  %>%
  mutate(treatment = gsub("Republican x ", "", var))

# democrat moderator
demhet_results <- comb_out_clean %>%
  filter(grepl("Democrat x", var)) %>%
  mutate(treatment = gsub("Democrat x ", "", var))

# independent moderator
indephet_results <- comb_out_clean %>%
  filter(grepl("Independent x", var)) %>%
  mutate(treatment = gsub("Independent x ", "", var))

# combine together
pidhet_results <- bind_rows("Republican" = rephet_results,
                            "Democrat" = demhet_results,
                            "Independent" = indephet_results,
                            .id = "Moderator")

pidhet_results <- pidhet_results %>%
  mutate(study_lab = case_when(study == "id17" ~ "Stokes and Warshaw (2017)",
                               study == "id35" ~ "Christenson and Goldfarb (2017)",
                               study == "id41" ~ "Singh and Swanson (2017)",
                               study == "id57" ~ "Schuldt and Enns (2017)",
                               study == "id71" ~ "Bolsen and Leeper (2014)",
                               study == "id73" ~ "Saunders (2017)",
                               study == "id74" ~ "Hardisty and Johnson (2010)",
                               study == "id83" ~ "Bolsen and Druckmann (2018)",
                               study == "id86" ~ "Schuldt, Konrath, and Schwarz (2011)",
                               study == "id115" ~ "DeGolia and Hiroyasu (2019)"))

pidhet_results <- pidhet_results %>%
  mutate(treatment_lab = str_replace(treatment, "treat: ", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization1b: ", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization2b: ", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization3b: ", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization4b: ", "")) %>%
  mutate(treatment_lab = str_replace(treatment_lab, "Randomization5b: ", ""))

# clean up
plot_pid <- pidhet_results %>%
  ggplot(aes(x = treatment_lab, y = `50%`)) +
  facet_grid(study ~ Moderator, scales = "free_y") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_pointrange(aes(ymin = `2.5%`, ymax = `97.5%`),
                  position = position_dodge(width = 1/2)) +
  # geom_point(aes(y = Mean), colour = "red") +
  ylab("Framing Effect") +
  xlab("Treatment") +
  coord_flip() + 
  theme_minimal() 

ggsave("./02_Paper/plot_pid.pdf", plot_pid, height = 1.7 * 7, width = 9)

# w/ OLS graphs

ols_comb_out %>%
  ggplot(aes(x = treatment_lab, y = estimate)) +
  facet_grid(study ~ Moderator, scales = "free_y") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 1/2)) +
  # geom_point(aes(y = Mean), colour = "red") +
  ylab("Framing Effect") +
  xlab("Treatment") +
  coord_flip() + 
  theme_minimal() 


pid_lasso_ols <- pidhet_results %>%
  ggplot(aes(x = treatment_lab, y = `50%`)) +
  facet_grid(study ~ Moderator, scales = "free_y") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_pointrange(aes(ymin = `2.5%`, ymax = `97.5%`),
                  position = position_dodge(width = 1/2)) +
  geom_pointrange(data = ols_comb_out,
                  aes(y = estimate, ymin = conf.low, ymax = conf.high,
                      colour = Moderator),
                  alpha = I(0.5)) +
  # geom_point(aes(y = Mean), colour = "red") +
  ylab("Framing Effect") +
  xlab("Treatment") +
  coord_flip() + 
  theme_minimal()

ggsave("./02_Paper/plot_pid_lasso_ols.pdf", pid_lasso_ols, height = 1.7 * 7, width = 9)

# alternative plot: scatterplot of OLS subgroup vs Lasso

pid_lasso_ols_comb <- ols_comb_out %>%
  dplyr::select(study, Moderator, outcome, treatment_lab, 
                estimate, conf.low, conf.high) %>%
  left_join(pidhet_results %>%
              dplyr::select(study, Moderator, outcome, treatment_lab, study_lab,
                            `50%`, `2.5%`, `97.5%`))

plot_scatter_pid_lasso_ols <- pid_lasso_ols_comb %>%
  ggplot(aes(x = estimate, y = `50%`,
             group = Moderator, colour = Moderator)) +
  facet_wrap( ~ Moderator, ncol = 2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_abline() +
  geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`),
                alpha = 0.2) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
                 alpha = 0.2) +
  geom_point() +
  geom_text(data = pid_lasso_ols_comb %>% filter(`50%` > 0.2),
            aes(label = study_lab, 
                x = estimate, y = `50%` + 0.3), size = 2) +
  coord_fixed(ratio = 1) +
  scale_colour_manual(values = c("Democrat" = "blue", "Independent" = "grey",
                                 "Republican" = "red")) +
  theme_minimal() +
  theme(legend.position = "none") +
  ylab("Sub-Group Effect (LASSOPlus)") +
  xlab("Sub-Group Effect (OLS)")

ggsave("./02_Paper/plot_scatter_pid_lasso_ols.pdf", plot_scatter_pid_lasso_ols, 
       height = 6, width = 1.7 * 6)

plot_scatter_pid_lasso_ols <- pid_lasso_ols_comb %>%
  ggplot(aes(x = estimate, y = `50%`,
             group = Moderator, colour = Moderator)) +
  facet_wrap( ~ Moderator, ncol = 1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_abline() +
  geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`),
                alpha = 0.2) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
                 alpha = 0.2) +
  geom_point() +
  geom_text(data = pid_lasso_ols_comb %>% filter(`50%` > 0.2),
            aes(label = study_lab, 
                x = estimate, y = `50%` + 0.3), size = 2) +
  coord_fixed(ratio = 1) +
  scale_colour_manual(values = c("Democrat" = "blue", "Independent" = "grey",
                                 "Republican" = "red")) +
  theme_minimal() +
  theme(legend.position = "none") +
  ylab("Sub-Group Effect (LASSOPlus)") +
  xlab("Sub-Group Effect (OLS)") + 
  coord_cartesian(ylim=c(-1, 1), xlim=c(-1, 1))

ggsave("./02_Paper/plot_scatter_pid_lasso_ols_rlim.pdf", plot_scatter_pid_lasso_ols, 
       height = 1.7 * 6, width =  6)


